Laplace distribution: laplace (double exponential)#
The Laplace distribution is a symmetric, continuous distribution with a sharp peak at its location and exponentially decaying tails. Compared to a Gaussian, it puts more mass near the center and more mass in the tails, which makes it a common choice for robust modeling.
In SciPy it appears as scipy.stats.laplace.
Learning goals#
understand what the Laplace distribution models and when it is useful
write down the PDF/CDF and connect them to sampling and likelihood
compute key moments (mean/variance/skewness/kurtosis) and entropy
derive the closed-form MLE (median + mean absolute deviation)
implement NumPy-only sampling and validate everything by Monte Carlo
Notebook roadmap#
Title & Classification
Intuition & Motivation
Formal Definition
Moments & Properties
Parameter Interpretation
Derivations
Sampling & Simulation (NumPy-only)
Visualization (PDF/CDF + Monte Carlo)
SciPy Integration
Statistical Use Cases
Pitfalls
Summary
import numpy as np
import plotly
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
import scipy
from scipy.stats import chi2, laplace, norm
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=4, suppress=True)
SEED = 7
rng = np.random.default_rng(SEED)
print("NumPy ", np.__version__)
print("SciPy ", scipy.__version__)
print("Plotly", plotly.__version__)
NumPy 1.26.2
SciPy 1.15.0
Plotly 6.5.2
1) Title & Classification#
Name:
laplaceType: continuous distribution
Support: x ∈ (-∞, ∞)
Parameter space: location μ ∈ ℝ and scale b > 0
We write:
The standard Laplace is (\mathrm{Laplace}(0, 1)).
2) Intuition & Motivation#
What it models#
The Laplace distribution is a natural model for real-valued errors that are:
centered around a typical value (μ)
more sharply peaked than a Gaussian
heavy-tailed relative to a Gaussian, but still with exponential tails
A clean way to remember it:
Normal: log-density is quadratic in ((x-\mu))
Laplace: log-density is linear in (|x-\mu|)
So Laplace errors correspond to an L1 loss.
Typical real-world use cases#
Robust regression / signal processing: modeling residuals with Laplace leads to least-absolute-deviation fitting.
Sparsity-inducing priors: the Laplace prior underpins the Bayesian view of Lasso (L1 regularization).
Differential privacy (Laplace mechanism): adding Laplace noise to numeric queries for privacy guarantees.
Noise with occasional big jumps: exponential tails can be a better match than Gaussian in some domains.
Relations to other distributions#
If (E_1, E_2) are i.i.d. exponential with mean (b), then (E_1 - E_2 \sim \mathrm{Laplace}(0, b)).
If (S\in{-1,+1}) is a fair sign and (E\sim\mathrm{Exp}(\text{mean}=b)), then (S,E\sim\mathrm{Laplace}(0, b)).
If (X\sim\mathrm{Laplace}(\mu, b)), then (|X-\mu|\sim\mathrm{Exp}(\text{mean}=b)).
The Laplace has heavier tails than a Gaussian but all moments exist; its MGF exists only for (|t|<1/b).
3) Formal Definition#
PDF#
For (X\sim\mathrm{Laplace}(\mu, b)):
CDF#
The CDF has a simple piecewise form:
Quantile function (inverse CDF)#
For (p\in(0,1)):
def _check_scale(b: float) -> float:
b = float(b)
if not np.isfinite(b) or b <= 0:
raise ValueError("`b` (scale) must be a positive, finite number.")
return b
def laplace_pdf(x, mu: float = 0.0, b: float = 1.0) -> np.ndarray:
b = _check_scale(b)
x = np.asarray(x, dtype=float)
return np.exp(-np.abs(x - mu) / b) / (2.0 * b)
def laplace_logpdf(x, mu: float = 0.0, b: float = 1.0) -> np.ndarray:
b = _check_scale(b)
x = np.asarray(x, dtype=float)
return -np.log(2.0 * b) - np.abs(x - mu) / b
def laplace_cdf(x, mu: float = 0.0, b: float = 1.0) -> np.ndarray:
b = _check_scale(b)
x = np.asarray(x, dtype=float)
z = (x - mu) / b
return np.where(x < mu, 0.5 * np.exp(z), 1.0 - 0.5 * np.exp(-z))
def laplace_ppf(p, mu: float = 0.0, b: float = 1.0, eps: float = 1e-12) -> np.ndarray:
b = _check_scale(b)
p = np.asarray(p, dtype=float)
if np.any((p <= 0) | (p >= 1)):
raise ValueError("p must be in (0, 1)")
p = np.clip(p, eps, 1.0 - eps)
left = mu + b * np.log(2.0 * p)
right = mu - b * (np.log(2.0) + np.log1p(-p)) # log(2(1-p)) stably
return np.where(p < 0.5, left, right)
# Quick cross-check vs SciPy
mu, b = 0.4, 1.7
x = np.linspace(-3, 3, 9)
p = np.array([0.1, 0.5, 0.9])
print("max |pdf - scipy|:", float(np.max(np.abs(laplace_pdf(x, mu, b) - laplace.pdf(x, loc=mu, scale=b)))))
print("max |cdf - scipy|:", float(np.max(np.abs(laplace_cdf(x, mu, b) - laplace.cdf(x, loc=mu, scale=b)))))
print("max |ppf - scipy|:", float(np.max(np.abs(laplace_ppf(p, mu, b) - laplace.ppf(p, loc=mu, scale=b)))))
max |pdf - scipy|: 0.0
max |cdf - scipy|: 0.0
max |ppf - scipy|: 0.0
4) Moments & Properties#
For (X\sim\mathrm{Laplace}(\mu, b)):
Mean: (\mathbb{E}[X]=\mu)
Variance: (\mathrm{Var}(X)=2b^2) (so (\mathrm{sd}(X)=\sqrt{2},b))
Skewness: 0 (symmetry)
Kurtosis: 6 (Pearson), so excess kurtosis = 3
MGF and characteristic function#
MGF exists only for (|t|<1/b):
Characteristic function (exists for all real (t)):
Entropy#
The differential entropy is:
def laplace_moments(mu: float = 0.0, b: float = 1.0):
b = _check_scale(b)
mean = float(mu)
var = float(2.0 * b * b)
skew = 0.0
kurt_excess = 3.0
entropy = float(1.0 + np.log(2.0 * b))
return mean, var, skew, kurt_excess, entropy
def sample_moments(x: np.ndarray):
x = np.asarray(x, dtype=float)
m = float(x.mean())
c = x - m
v = float(np.mean(c**2))
skew = float(np.mean(c**3) / (v ** 1.5))
kurt_excess = float(np.mean(c**4) / (v**2) - 3.0)
return m, v, skew, kurt_excess
mu, b = 1.5, 0.8
mean_f, var_f, skew_f, kurt_f, ent_f = laplace_moments(mu, b)
mean_s, var_s, skew_s, kurt_s = laplace.stats(loc=mu, scale=b, moments="mvsk")
ent_s = laplace.entropy(loc=mu, scale=b)
print("theory (mean, var, skew, kurt_excess, entropy):", (mean_f, var_f, skew_f, kurt_f, ent_f))
print("scipy (mean, var, skew, kurt_excess, entropy):", (float(mean_s), float(var_s), float(skew_s), float(kurt_s), float(ent_s)))
x = laplace.rvs(loc=mu, scale=b, size=300_000, random_state=rng)
mean_mc, var_mc, skew_mc, kurt_mc = sample_moments(x)
print("monte (mean, var, skew, kurt_excess): ", (mean_mc, var_mc, skew_mc, kurt_mc))
theory (mean, var, skew, kurt_excess, entropy): (1.5, 1.2800000000000002, 0.0, 3.0, 1.4700036292457357)
scipy (mean, var, skew, kurt_excess, entropy): (1.5, 1.2800000000000002, 0.0, 3.0, 1.4700036292457357)
monte (mean, var, skew, kurt_excess): (1.501420070672676, 1.2685335280033816, -0.0003665014939261965, 2.9428857779033946)
def laplace_cf(t, mu: float = 0.0, b: float = 1.0) -> np.ndarray:
b = _check_scale(b)
t = np.asarray(t, dtype=float)
return np.exp(1j * mu * t) / (1.0 + (b * t) ** 2)
mu, b = 0.0, 1.0
t = np.linspace(-15, 15, 2000)
phi = laplace_cf(t, mu=mu, b=b)
fig = make_subplots(rows=1, cols=2, subplot_titles=("Re φ(t)", "Im φ(t)"))
fig.add_trace(go.Scatter(x=t, y=np.real(phi), mode="lines"), row=1, col=1)
fig.add_trace(go.Scatter(x=t, y=np.imag(phi), mode="lines"), row=1, col=2)
fig.update_xaxes(title_text="t", row=1, col=1)
fig.update_xaxes(title_text="t", row=1, col=2)
fig.update_layout(width=950, height=350, showlegend=False, title="Characteristic function of Laplace(0,1)")
fig.show()
5) Parameter Interpretation#
μ (location) shifts the distribution left/right. For Laplace it is simultaneously the mean, median, and mode.
b (scale) controls both the peak height and tail thickness:
(f(\mu)=1/(2b)) so larger b means a lower peak.
(\mathrm{sd}(X)=\sqrt{2},b) so b is proportional to standard deviation.
Below we visualize how ((\mu,b)) changes the PDF and CDF.
x = np.linspace(-10, 10, 3000)
params = [
(0.0, 0.5),
(0.0, 1.0),
(0.0, 2.0),
(2.0, 1.0),
]
fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))
for mu, b in params:
fig.add_trace(
go.Scatter(x=x, y=laplace_pdf(x, mu=mu, b=b), mode="lines", name=f"μ={mu}, b={b}"),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(x=x, y=laplace_cdf(x, mu=mu, b=b), mode="lines", name=f"μ={mu}, b={b}"),
row=1,
col=2,
)
fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(width=1000, height=420)
fig.show()
6) Derivations#
Expectation#
Let (Y=X-\mu). Then (Y) has symmetric density (f_Y(y) = (1/(2b))\exp(-|y|/b)).
Because (y f_Y(y)) is an odd function,
Variance#
Using symmetry again,
The integral is a Gamma-type integral: (\int_0^{\infty} y^2 e^{-y/b}dy = 2!,b^3 = 2b^3), hence
Likelihood (and the MLE)#
For i.i.d. data (x_1,\dots,x_n) from (\mathrm{Laplace}(\mu,b)), the likelihood is
So the log-likelihood is
Key consequence: for fixed (b), maximizing (\ell) over (\mu) is equivalent to minimizing (\sum_i |x_i-\mu|), which is minimized by any median of the sample.
With (\hat\mu) chosen as a sample median, differentiating (\ell(\hat\mu,b)) with respect to (b) gives the (closed-form) MLE:
def laplace_loglik(data: np.ndarray, mu: float, b: float) -> float:
return float(np.sum(laplace_logpdf(data, mu=mu, b=b)))
def laplace_mle_closed_form(data: np.ndarray, b_floor: float = 1e-12) -> tuple[float, float]:
"""Closed-form MLE for Laplace(μ,b): μ=median, b=mean absolute deviation from μ.
Notes
-----
- If all observations are identical, the likelihood increases as b -> 0+.
We return a small positive floor to stay inside the parameter space.
"""
x = np.asarray(data, dtype=float)
if x.ndim != 1:
raise ValueError("data must be 1D")
mu_hat = float(np.median(x))
b_hat = float(np.mean(np.abs(x - mu_hat)))
b_hat = float(max(b_hat, b_floor))
return mu_hat, b_hat
mu_true, b_true = 1.0, 0.7
data = laplace.rvs(loc=mu_true, scale=b_true, size=2000, random_state=rng)
mu_hat, b_hat = laplace_mle_closed_form(data)
mu_hat_sp, b_hat_sp = laplace.fit(data)
print("true (μ, b) =", (mu_true, b_true))
print("closed-form MLE =", (mu_hat, b_hat))
print("scipy laplace.fit =", (float(mu_hat_sp), float(b_hat_sp)))
print("loglik at MLE =", laplace_loglik(data, mu_hat, b_hat))
true (μ, b) = (1.0, 0.7)
closed-form MLE = (1.000653645441159, 0.716006489989578)
scipy laplace.fit = (1.000653645441159, 0.716006489989578)
loglik at MLE = -2718.1622654572566
7) Sampling & Simulation (NumPy-only)#
Two convenient sampling views:
Inverse CDF: sample (U\sim\mathrm{Uniform}(0,1)) and return (F^{-1}(U)).
Difference of exponentials: if (E_1,E_2\overset{iid}{\sim}\mathrm{Exp}(\text{mean}=b)), then (\mu + (E_1-E_2)\sim\mathrm{Laplace}(\mu,b)).
We’ll implement both with NumPy and sanity-check them.
def laplace_rvs_inverse(
rng: np.random.Generator,
size: int,
mu: float = 0.0,
b: float = 1.0,
eps: float = 1e-12,
) -> np.ndarray:
b = _check_scale(b)
u = rng.random(size)
u = np.clip(u, eps, 1.0 - eps)
left = mu + b * np.log(2.0 * u)
right = mu - b * (np.log(2.0) + np.log1p(-u))
return np.where(u < 0.5, left, right)
def laplace_rvs_exp_difference(
rng: np.random.Generator,
size: int,
mu: float = 0.0,
b: float = 1.0,
) -> np.ndarray:
b = _check_scale(b)
e1 = rng.exponential(scale=b, size=size)
e2 = rng.exponential(scale=b, size=size)
return mu + (e1 - e2)
mu, b = -0.5, 1.2
n = 200_000
x_inv = laplace_rvs_inverse(rng, size=n, mu=mu, b=b)
x_diff = laplace_rvs_exp_difference(rng, size=n, mu=mu, b=b)
q = [0.05, 0.25, 0.5, 0.75, 0.95]
print("theory quantiles:", np.round(laplace.ppf(q, loc=mu, scale=b), 4))
print("inv quantiles:", np.round(np.quantile(x_inv, q), 4))
print("diff quantiles:", np.round(np.quantile(x_diff, q), 4))
theory quantiles: [-3.2631 -1.3318 -0.5 0.3318 2.2631]
inv quantiles: [-3.2641 -1.3348 -0.5003 0.3281 2.2367]
diff quantiles: [-3.2562 -1.3277 -0.4983 0.3328 2.2574]
8) Visualization (PDF/CDF + Monte Carlo)#
We’ll plot the theoretical PDF/CDF and compare them to Monte Carlo samples.
mu, b = 0.0, 1.0
x = np.linspace(-8, 8, 4000)
fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))
fig.add_trace(go.Scatter(x=x, y=laplace_pdf(x, mu=mu, b=b), mode="lines", name="pdf"), row=1, col=1)
fig.add_trace(go.Scatter(x=x, y=laplace_cdf(x, mu=mu, b=b), mode="lines", name="cdf"), row=1, col=2)
fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(width=950, height=380, showlegend=False, title="Laplace(0,1): PDF and CDF")
fig.show()
# Monte Carlo histogram vs theory
n = 120_000
samples = laplace_rvs_inverse(rng, size=n, mu=mu, b=b)
fig = go.Figure()
fig.add_trace(
go.Histogram(
x=samples,
histnorm="probability density",
nbinsx=120,
name="Monte Carlo",
opacity=0.6,
)
)
fig.add_trace(go.Scatter(x=x, y=laplace_pdf(x, mu=mu, b=b), mode="lines", name="Theory"))
fig.update_layout(
title="Laplace(0,1): histogram vs PDF",
xaxis_title="x",
yaxis_title="density",
width=900,
height=420,
barmode="overlay",
)
fig.show()
# Empirical CDF vs theory
xs = np.sort(samples)
ecdf = np.arange(1, xs.size + 1) / xs.size
fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=ecdf, mode="lines", name="Empirical"))
fig.add_trace(go.Scatter(x=x, y=laplace_cdf(x, mu=mu, b=b), mode="lines", name="Theory"))
fig.update_layout(title="Laplace(0,1): empirical CDF", xaxis_title="x", yaxis_title="F(x)", width=900, height=420)
fig.show()
# Q-Q plot (sample quantiles vs theoretical quantiles)
p = np.linspace(0.01, 0.99, 200)
q_theory = laplace_ppf(p, mu=mu, b=b)
q_emp = np.quantile(samples, p)
fig = go.Figure()
fig.add_trace(go.Scatter(x=q_theory, y=q_emp, mode="markers", name="quantiles"))
fig.add_trace(go.Scatter(x=q_theory, y=q_theory, mode="lines", name="y=x"))
fig.update_layout(title="Laplace(0,1): Q-Q plot", xaxis_title="theoretical quantile", yaxis_title="sample quantile", width=700, height=520)
fig.show()
9) SciPy Integration#
SciPy’s laplace uses the same parameterization:
loc= μscale= b
The most common methods:
laplace.pdf(x, loc, scale)/laplace.logpdf(...)laplace.cdf(x, loc, scale)/laplace.ppf(p, loc, scale)laplace.rvs(loc, scale, size, random_state=...)laplace.fit(data)(MLE)
mu, b = 0.7, 1.2
x = np.array([-2.0, 0.0, 1.0, 3.0])
print("pdf:", laplace.pdf(x, loc=mu, scale=b))
print("cdf:", laplace.cdf(x, loc=mu, scale=b))
print("rvs:", laplace.rvs(loc=mu, scale=b, size=5, random_state=rng))
data = laplace.rvs(loc=mu, scale=b, size=3000, random_state=rng)
mu_hat_sp, b_hat_sp = laplace.fit(data)
mu_hat_cf, b_hat_cf = laplace_mle_closed_form(data)
print("fit (scipy) :", (float(mu_hat_sp), float(b_hat_sp)))
print("fit (closed-form):", (mu_hat_cf, b_hat_cf))
pdf: [0.0439 0.2325 0.3245 0.0613]
cdf: [0.0527 0.279 0.6106 0.9265]
rvs: [1.3046 0.8555 4.5273 1.1789 0.7094]
fit (scipy) : (0.7282358708034962, 1.1857197356240277)
fit (closed-form): (0.7282358708034962, 1.1857197356240277)
10) Statistical Use Cases#
Hypothesis testing#
The Laplace likelihood leads to simple likelihood-ratio tests (LRTs). We’ll demonstrate an LRT for the location parameter:
(H_0: \mu = \mu_0)
(H_1: \mu) free
We keep (b) unknown in both models (estimated by MLE). Asymptotically, the LRT statistic is (\chi^2) with 1 degree of freedom.
Bayesian modeling#
A Laplace prior (p(\theta) \propto \exp(-|\theta|/b)) induces sparsity/shrinkage. Under a Gaussian likelihood, the MAP estimator becomes soft-thresholding (the 1D analog of Lasso).
Generative modeling#
Laplace noise is used as a primitive generator for privacy-preserving releases (Laplace mechanism), and also as a building block in mixture models and robust generative pipelines.
def laplace_lrt_mu(data: np.ndarray, mu0: float) -> tuple[float, float]:
"""Likelihood-ratio test for H0: μ=mu0 vs H1: μ free (b unknown in both)."""
x = np.asarray(data, dtype=float)
mu_hat, b_hat = laplace_mle_closed_form(x)
b0_hat = float(np.mean(np.abs(x - mu0)))
b0_hat = float(max(b0_hat, 1e-12))
ll_hat = laplace_loglik(x, mu=mu_hat, b=b_hat)
ll_0 = laplace_loglik(x, mu=mu0, b=b0_hat)
LR = 2.0 * (ll_hat - ll_0)
p_value = float(chi2.sf(LR, df=1))
return float(LR), p_value
mu0 = 0.0
b_true = 1.0
n = 120
# Under H0
data0 = laplace_rvs_inverse(rng, size=n, mu=mu0, b=b_true)
LR0, p0 = laplace_lrt_mu(data0, mu0=mu0)
print(f"H0 sample: LR={LR0:.3f}, p={p0:.3f}")
# Under H1
mu_true = 0.7
data1 = laplace_rvs_inverse(rng, size=n, mu=mu_true, b=b_true)
LR1, p1 = laplace_lrt_mu(data1, mu0=mu0)
print(f"H1 sample: LR={LR1:.3f}, p={p1:.3e}")
# Quick type-I check (asymptotic calibration)
m = 1500
pvals = np.empty(m)
for i in range(m):
d = laplace_rvs_inverse(rng, size=n, mu=mu0, b=b_true)
_, pvals[i] = laplace_lrt_mu(d, mu0=mu0)
print("approx P(p<0.05) under H0:", float(np.mean(pvals < 0.05)))
H0 sample: LR=0.151, p=0.698
H1 sample: LR=60.434, p=7.608e-15
approx P(p<0.05) under H0: 0.056666666666666664
def soft_threshold(y: np.ndarray, tau: float) -> np.ndarray:
y = np.asarray(y, dtype=float)
return np.sign(y) * np.maximum(np.abs(y) - tau, 0.0)
# Bayesian demo: Gaussian likelihood with Laplace prior
# y | theta ~ N(theta, sigma^2), theta ~ Laplace(0, b)
sigma = 1.0
b_prior = 0.7
# MAP has a closed-form soft-threshold in this 1D model
tau = sigma**2 / b_prior
ys = np.linspace(-4, 4, 401)
theta_map = soft_threshold(ys, tau=tau)
fig = go.Figure()
fig.add_trace(go.Scatter(x=ys, y=theta_map, mode="lines", name="MAP( y )"))
fig.add_trace(go.Scatter(x=ys, y=ys, mode="lines", name="MLE=y", line=dict(dash="dash")))
fig.update_layout(
title=f"Laplace prior shrinkage (b={b_prior}, sigma={sigma})",
xaxis_title="observation y",
yaxis_title="MAP estimate of theta",
width=850,
height=420,
)
fig.show()
# Visualize the posterior for a single observation
y0 = 1.2
theta_grid = np.linspace(-5, 5, 2001)
log_like = norm.logpdf(y0, loc=theta_grid, scale=sigma)
log_prior = laplace_logpdf(theta_grid, mu=0.0, b=b_prior)
log_post = log_like + log_prior
log_post -= log_post.max()
post_unnorm = np.exp(log_post)
post = post_unnorm / np.trapz(post_unnorm, theta_grid)
theta_map_grid = float(theta_grid[np.argmax(post)])
theta_map_closed = float(soft_threshold(y0, tau=tau))
fig = go.Figure()
fig.add_trace(go.Scatter(x=theta_grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=theta_map_grid, line_dash="dash", line_color="black", annotation_text="MAP")
fig.update_layout(
title=f"Posterior p(theta|y={y0}) with Laplace prior",
xaxis_title="theta",
yaxis_title="density",
width=900,
height=420,
)
fig.show()
print("MAP (grid) :", theta_map_grid)
print("MAP (closed):", theta_map_closed)
MAP (grid) : 0.0
MAP (closed): 0.0
# Generative / differential privacy-style demo: Laplace mechanism for a count query
# For a count, sensitivity Δf = 1. Laplace mechanism adds noise ~ Laplace(0, Δf/ε) = Laplace(0, 1/ε).
true_count = 250
eps_values = [0.25, 0.5, 1.0, 2.0] # larger ε => less noise
n = 60_000
fig = go.Figure()
for eps in eps_values:
b_noise = 1.0 / eps
noise = laplace_rvs_inverse(rng, size=n, mu=0.0, b=b_noise)
released = true_count + noise
fig.add_trace(
go.Histogram(
x=released,
histnorm="probability density",
nbinsx=140,
opacity=0.35,
name=f"ε={eps} (b={b_noise:.2f})",
)
)
fig.update_layout(
title="Laplace mechanism: noisy releases of a count",
xaxis_title="released value",
yaxis_title="density",
barmode="overlay",
width=950,
height=450,
)
fig.show()
# Tail bound: P(|Noise| > t) = exp(-t/b) = exp(-ε t)
alpha = 0.05
for eps in eps_values:
b_noise = 1.0 / eps
t95 = b_noise * np.log(1.0 / alpha)
print(f"ε={eps:>4}: P(|err|>t)={alpha} at t={t95:.3f}")
ε=0.25: P(|err|>t)=0.05 at t=11.983
ε= 0.5: P(|err|>t)=0.05 at t=5.991
ε= 1.0: P(|err|>t)=0.05 at t=2.996
ε= 2.0: P(|err|>t)=0.05 at t=1.498
11) Pitfalls#
Invalid parameters: the scale must satisfy b > 0. Many routines will return NaNs or raise errors if b ≤ 0.
Scale vs standard deviation: for Laplace, (\mathrm{sd}=\sqrt{2},b). If you want a target variance (\sigma^2), set (b=\sigma/\sqrt{2}).
MGF domain: (M_X(t)) only exists for (|t|<1/b). (The characteristic function exists for all t.)
Sampling numerics: inverse-CDF sampling involves (\log U). Clip
Uaway from 0 and 1 to avoidlog(0).Degenerate samples: if all observations are identical, the likelihood increases as b → 0+. Any practical MLE needs a small floor.
Fitting and medians: for even n, the minimizing set of medians is an interval.
np.medianreturns the midpoint, which is still optimal.
12) Summary#
laplaceis a continuous distribution on ((-∞,∞)) with location μ and scale b>0.PDF: (f(x)=\frac{1}{2b}\exp(-|x-\mu|/b)); CDF is piecewise exponential.
Mean = μ, variance = (2b^2), skewness = 0, excess kurtosis = 3, entropy = (1+\log(2b)).
MLE: (\hat\mu) is a sample median and (\hat b) is the mean absolute deviation from (\hat\mu).
Sampling (NumPy-only): inverse CDF or difference of exponentials.
Common uses: robust modeling (L1 loss), sparse priors (Bayesian Lasso), and privacy noise (Laplace mechanism).